﻿#include "discrete_face_spline.hxx"
#include "discrete_face_loop.hxx"
#include "idx_mesh_utils.hxx"
#include "KDTree.hpp"
// #define DEBUG
#ifdef DEBUG
#    define TIMERSTART(tag) auto tag##_start = std::chrono::high_resolution_clock::now();
#else
#    define TIMERSTART(tag)
#endif  // DEBUG

#ifdef DEBUG
#    define TIMEREND(tag)                                           \
        auto tag##_end = std::chrono::high_resolution_clock::now(); \
        printf("|| %-40s|| time costs %f ms\n", #tag, std::chrono::duration<double, std::milli>(tag##_end - tag##_start).count());
#else
#    define TIMEREND(tag)
#endif  // DEBUG

QuadTree::QuadTree(FACE* f, const Box box, double nt, double st, std::vector<SPApar_pos> vertex_set = std::vector<SPApar_pos>(0)) : face_(f), mBox_(box), distance_tolerance_(st), angle_tolerance_(nt), mRoot_(std::make_unique<Node>()) {
    // 获取主参数域范围
    const surface& surf = f->geometry()->equation();
    // 根据不同情况分割
    // 如果角点个数为1，则最开始需要四叉划分
    // SPAposition corners_pos_3d[4] = {surf.eval_position(mBox_.getBotLeft()), surf.eval_position(mBox_.getTopLeft()), surf.eval_position(mBox_.getTopRight()), surf.eval_position(mBox_.getBotRight())};
    // int corners_num = 4;
    if (surf.closed_u() && surf.closed_v()) {
        need_quad_first_ = true;
    }
    else if (surf.closed_u() && surf.singular_v(surf.param_range_v().start_pt()) && surf.singular_v(surf.param_range_v().end_pt())) {
        need_quad_first_ = true;
    }
    else if (surf.closed_v() && surf.singular_u(surf.param_range_u().start_pt()) && surf.singular_u(surf.param_range_u().end_pt())) {
        need_quad_first_ = true;
    }
    // 通过二叉划分构造树，并保存所有的叶子节点边界框
    BuildTreeByBin(mBox_, 1);
    // 从所有的叶子节点边界框中生成初始采样点和边
    GenaratePointsAndEdges();
}
QuadTree::~QuadTree() {
    for (auto x : v_iso_lines) {
        ACIS_DELETE x.second;
    }
    for (auto x : u_iso_lines) {
        ACIS_DELETE x.second;
    }
}

void QuadTree::GenaratePointsAndEdges() {
    using Edge = std::pair<int, int>;
    // 定义四个角点坐标的索引
    struct CornerIndexes {
        int bot_left;
        int bot_right;
        int top_left;
        int top_right;
    };

    std::unordered_map<point, int> point2index_map;  // 点到索引的映射，相当于维护一个去重点集
    int point_count = 0;                             // 点的计数
    std::vector<CornerIndexes> index_of_corners;     // 四个角点的索引
    // 获取或创建点的索引
    auto GetOrCreatePointIndex = [&](const SPApar_pos& spa_pos) {
        const point p = { spa_pos.u, spa_pos.v };                 // 将参数域坐标转换为点
        if (point2index_map.find(p) == point2index_map.end()) {  // 如果点不在映射中，说明是新的点
            point2index_map[p] = point_count++;                 // 将点添加到映射中，并分配新的索引
            points_.push_back(p);                               // 将点添加到点的容器中
            return point_count - 1;                             // 返回新分配的索引
        }
        return point2index_map[p];  // 如果点已经存在，则直接返回已有的索引
        };

    index_of_corners.reserve(leaf_rect_.size());
    // 遍历所有叶节点，生成矩形的四个顶点的索引
    for (int i = 0; i < leaf_rect_.size(); i++) {
        int bot_left_idx = GetOrCreatePointIndex(leaf_rect_[i].getBotLeft());                     // 左下角
        int bot_right_idx = GetOrCreatePointIndex(leaf_rect_[i].getBotRight());                   // 右下角
        int top_right_idx = GetOrCreatePointIndex(leaf_rect_[i].getTopRight());                   // 右上角
        int top_left_idx = GetOrCreatePointIndex(leaf_rect_[i].getTopLeft());                     // 左上角
        index_of_corners.emplace_back(bot_left_idx, bot_right_idx, top_left_idx, top_right_idx);  // 将矩形的四个顶点索引添加到容器中
    }

    // 并查集合并短边
    std::vector<int> box_indices(leaf_rect_.size());       // 创建边界框的索引，用于边界框的排序，我们不直接对边界框进行排序，而是通过索引间接排序
    std::iota(box_indices.begin(), box_indices.end(), 0);  // 初始化边界框索引，从0到n-1
    std::vector<int> parent(point_count);                  // 并查集的父节点数组
    std::function<int(int)> find;                          // 并查集的查找函数
    find = [&](int x) -> int {
        if (x != parent[x]) {
            parent[x] = find(parent[x]);  // 路径压缩
        }
        return parent[x];
        };

    // 此lambda用于处理边界框的边，根据边界框的宽度或者高度排序，长度升序地处理边，还可以过滤长边
    auto process_edges = [&](auto comparator, auto get_edge1, auto get_edge2, edge_vector& edges) {
        std::sort(box_indices.begin(), box_indices.end(), comparator);  // 按照比较器排序边界框索引
        std::iota(parent.begin(), parent.end(), 0);                     // 初始化并查集
        for (int i = 0; i < box_indices.size(); i++) {
            const std::array<Edge, 2>& edges_to_check = { get_edge1(box_indices[i]), get_edge2(box_indices[i]) };  // 获取需要检查的边
            for (const auto& [idx_a, idx_b] : edges_to_check) {
                int root_a = find(idx_a);
                int root_b = find(idx_b);
                if (root_a != root_b) {                 // 如果根不同，说明是不同的集合，说明该边不会由其他更短的边组成
                    edges.emplace_back(idx_a, idx_b);  // 将边添加到边的容器中
                    parent[root_a] = root_b;           // 合并集合
                }
            }
        }
        };

    // 处理水平边
    process_edges([&](int a, int b) { return leaf_rect_[a].getWidth() < leaf_rect_[b].getWidth(); }, [&](int rect_id) { return std::make_pair(index_of_corners[rect_id].bot_left, index_of_corners[rect_id].bot_right); },
        [&](int rect_id) { return std::make_pair(index_of_corners[rect_id].top_left, index_of_corners[rect_id].top_right); }, h_edges_);

    // 处理垂直边
    process_edges([&](int a, int b) { return leaf_rect_[a].getHight() < leaf_rect_[b].getHight(); }, [&](int rect_id) { return std::make_pair(index_of_corners[rect_id].bot_left, index_of_corners[rect_id].top_left); },
        [&](int rect_id) { return std::make_pair(index_of_corners[rect_id].top_right, index_of_corners[rect_id].bot_right); }, v_edges_);
}

bool QuadTree::isLeaf(const Node* node) const {
    return !static_cast<bool>(node->children[0]);
}
Box QuadTree::GetRootBox() const {
    return mBox_;
}
void QuadTree::BuildTreeByBin(const Box& box, int dep) {
    // 判断终止条件
    if ((dep > 1 && CheckFitAndJudgeDiagonal(face_, box, angle_tolerance_, distance_tolerance_))) {  // Avoid error recursion with magic number on tapered tip
        leaf_rect_.push_back(box);                                                                  // 叶子节点，存下边界框
        return;
    }
    const surface& spl = face_->geometry()->equation();
    // 构造U和V向中线曲线，用于近似表示曲面曲率
    if (u_iso_lines.find(box.getCenter().v) == u_iso_lines.end()) {
        u_iso_lines[box.getCenter().v] = spl.u_param_line(box.getCenter().v);
    }
    if (v_iso_lines.find(box.getCenter().u) == v_iso_lines.end()) {
        v_iso_lines[box.getCenter().u] = spl.v_param_line(box.getCenter().u);
    }
    // 获取曲线的长度
    auto u_line = u_iso_lines[box.getCenter().v];
    double u_len = u_line->length(box.getLeft(), box.getRight());
    auto v_line = v_iso_lines[box.getCenter().u];
    double v_len = v_line->length(box.getBot(), box.getTop());
    // 获取两条曲线的两个端点之间的三维欧几里得距离
    double u_distance_3d = (spl.eval_position(box.getLeftMid()) - spl.eval_position(box.getRightMid())).len();
    double v_distance_3d = (spl.eval_position(box.getBotMid()) - spl.eval_position(box.getTopMid())).len();

    bool is_two_straight_line = is_zero_mch(u_len - u_distance_3d) && is_zero_mch(v_len - v_distance_3d);

    // 近似的曲率比率计算
    double ur = is_zero_mch(u_distance_3d) ? 0 : (u_len - u_distance_3d) / u_distance_3d;
    double vr = is_zero_mch(v_distance_3d) ? 0 : (v_len - v_distance_3d) / v_distance_3d;

    if (need_quad_first_ || (is_zero_mch(v_distance_3d) && is_zero_mch(u_distance_3d)) || is_two_straight_line || (!is_zero_mch(ur) && !is_zero_mch(vr) && is_equal(ur, vr))) {
        if (need_quad_first_) need_quad_first_ = false;
        BuildTreeByBin(Box{ box.getBotLeft(), box.getCenter() }, dep + 1);
        BuildTreeByBin(Box{ box.getBotMid(), box.getRightMid() }, dep + 1);
        BuildTreeByBin(Box{ box.getCenter(), box.getTopRight() }, dep + 1);
        BuildTreeByBin(Box{ box.getLeftMid(), box.getTopMid() }, dep + 1);
        return;
    }
    // 若u向为0 或者 u向比率大于v向比率，则按照u向划分，否则v向划分
    if (is_zero_mch(u_distance_3d) || (ur > vr && !is_zero_mch(v_distance_3d))) {
        BuildTreeByBin(Box{ box.getBotLeft(), box.getTopMid() }, dep + 1);
        BuildTreeByBin(Box{ box.getBotMid(), box.getTopRight() }, dep + 1);
    }
    else {
        BuildTreeByBin(Box{ box.getBotLeft(), box.getRightMid() }, dep + 1);
        BuildTreeByBin(Box{ box.getLeftMid(), box.getTopRight() }, dep + 1);
    }
    return;
}
void spline_faceter::init() {
    spl = (spline*)&f->geometry()->equation_for_update();
    approx_tol_ = std::max(spl->fitol(), approx_tol_);
    SPApar_box face_par_box = spl->param_range();
    SPAinterval u_interval = spl->param_range_u();
    SPAinterval v_interval = spl->param_range_v();
    closed_u = spl->closed_u();
    closed_v = spl->closed_v();
    periodic_u = spl->periodic_u();
    periodic_v = spl->periodic_v();
    u_le = u_interval.start_pt();
    u_ri = u_interval.end_pt();
    v_le = v_interval.start_pt();
    v_ri = v_interval.end_pt();
    singular_vle = spl->singular_u(v_le);
    singular_vri = spl->singular_u(v_ri);
    singular_ule = spl->singular_u(u_le);
    singular_uri = spl->singular_u(u_ri);
    u_tolerance = spl->param_range_u().length() / 2;
    v_tolerance = spl->param_range_v().length() / 2;
    if (spl->closed_u()) {
        if (spl->singular_u(v_le)) {
            singular_vle = true;
        }
        else {  // B_7  B14的样例需要通过下列判断才能确定奇点情况
            curve* vle = spl->v_param_line(v_le);
            SPAinterval interval = vle->param_range();
            if (is_equal(vle->length(interval.start_pt(), interval.end_pt()), 0.0)) {
                singular_vle = true;
            }
        }
        if (spl->singular_u(v_ri)) {
            singular_vri = true;
        }
        else {
            curve* vri = spl->v_param_line(v_ri);
            SPAinterval interval = vri->param_range();
            if (is_equal(vri->length(interval.start_pt(), interval.end_pt()), 0.0)) {
                singular_vri = true;
            }
        }
    }
    if (spl->closed_v()) {
        if (spl->singular_u(u_le)) {
            singular_ule = true;
        }
        else {  // B_7  B14的样例需要通过下列判断才能确定奇点情况
            curve* ule = spl->v_param_line(u_le);
            SPAinterval interval = ule->param_range();
            if (is_equal(ule->length(interval.start_pt(), interval.end_pt()), 0.0)) {
                singular_ule = true;
            }
        }
        if (spl->singular_u(u_ri)) {
            singular_uri = true;
        }
        else {
            curve* uri = spl->v_param_line(u_ri);
            SPAinterval interval = uri->param_range();
            if (is_equal(uri->length(interval.start_pt(), interval.end_pt()), 0.0)) {
                singular_uri = true;
            }
        }
    }
}

void spline_faceter::peripheryProcess() {
}

bool spline_faceter::USeperationProcess() {
    return false;  // 返回false
}

bool spline_faceter::VSeperationProcess() {
    return false;  // 返回false
}

void spline_faceter::holeProcess() {
    // 检查是否有孔洞Loop需要处理
    if (hole_loops.size()) {
        for (int i = 0; i < hole_loops.size(); i++) {
            if (hole_loops[i].size() < 2) continue;
            point_vector loop = std::move(hole_loops[i]);
            size_t start = points.size();
            for (int p = 0; p < loop.size(); p++) {
                points.push_back(loop[p]);
                edges.push_back(CDT::Edge(points.size() - 1, points.size()));
            }
            edges.pop_back();
            edges.push_back(CDT::Edge(points.size() - 1, start));
        }
    }
}

void spline_faceter::attachMesh() {
    SPAinterval u_interval = f->geometry()->equation().param_range_u();
    SPAinterval v_interval = f->geometry()->equation().param_range_v();
    // 初始化叶子节点的边界框
    Box box = Box{
      {u_interval.start_pt(), v_interval.start_pt()},
      {u_interval.end_pt(),   v_interval.end_pt()  }
    };
    QuadTree tree(f, box, nt, st);
    attachMeshByQuad(tree);
}

DeBugOutputUVGraph debug_graph;
void spline_faceter::ProcessFaceAndLoopsIntersection(QuadTree& tree, point_vector& loops_points, edge_vector& loops_edges, point_vector& res_points, edge_vector& res_edges, CDT::EdgeUSet& markedEdge) {
    const int NOTFOUND = -1;
    point_vector face_points = std::move(tree.points_);                                                                            // 面采样获得的点
    edge_vector v_edges = std::move(tree.v_edges_);                                                                                // 面采样得到的垂直边
    edge_vector h_edges = std::move(tree.h_edges_);                                                                                // 面采样得到的水平边
    auto v_intersection_results = CalculatePolygonEdgeIntersections(face_points, v_edges, loops_edges, loops_points, GetX, GetY);  // 计算垂直边与Loop的交点
    auto h_intersection_results = CalculatePolygonEdgeIntersections(face_points, h_edges, loops_edges, loops_points, GetY, GetX);  // 计算水平边与Loop的交点

    point_vector new_points;                                                                // 新点集
    std::vector<int> new_map(face_points.size(), NOTFOUND);                                 // 点索引映射，初始化为NOTFOUND
    auto res = JudgePointsInPolygon2D(face_points, loops_edges, loops_points, GetX, GetY);  // 判断点是否在多边形内部
    auto res_ = JudgePointsInPolygon2D(face_points, loops_edges, loops_points, GetY, GetX);
    // 返回的两个res是不在面内的点的集合，故做一个并集
    res.insert(res_.begin(), res_.end());
    for (int i = 0; i < face_points.size(); i++) {  // 遍历面的点
        if (res.find(i) == res.end()) {             // 如果点不在集合中，说明是面内点
            new_points.push_back(face_points[i]);  // 添加到新点集
            new_map[i] = new_points.size() - 1;    // 更新索引映射
        }
    }
    // 以下是使用KDtree查找最近点对
    pointVec tmp;
    for (auto& x : loops_points) {
        auto pos_3d = spl->eval_position({ x.x, x.y });
        tmp.push_back({ pos_3d.x(), pos_3d.y(), pos_3d.z() });
    }
    KDTree kdtree = KDTree(tmp);  // 创建KD树
    // 以下用于查找最近点,三维欧几里得距离第一关键词，二维参数域欧几里得距离第二关键词
    auto find_near_idx_from_kdtree = [&](point temp) -> int {
        SPAposition pos_3d = spl->eval_position({ temp.x, temp.y });  // 获取交点的三维坐标
        indexArr near_idxs = kdtree.nearest_indices({ pos_3d.x(), pos_3d.y(), pos_3d.z() }, 4);
        int near_idx = near_idxs.front();
        double min_len_sq = (spl->eval_position({ res_points[near_idx].x, res_points[near_idx].y }) - pos_3d).len_sq();
        for (int id : near_idxs) {
            double len_sq = (spl->eval_position({ res_points[id].x, res_points[id].y }) - pos_3d).len_sq();
            if (len_sq < min_len_sq) {
                min_len_sq = len_sq;
                near_idx = id;
            }
        }
        for (auto id : near_idxs) {  // 使用的是CDT的point类型，只能手写二维欧几里得距离函数
            if (ba.SamePointJudgeBy3D({ res_points[near_idx].x, res_points[near_idx].y }, { res_points[id].x, res_points[id].y }) &&
                (temp.x - res_points[near_idx].x) * (temp.x - res_points[near_idx].x) + (temp.y - res_points[near_idx].y) * (temp.y - res_points[near_idx].y) >
                (temp.x - res_points[id].x) * (temp.x - res_points[id].x) + (temp.y - res_points[id].y) * (temp.y - res_points[id].y)) {
                near_idx = id;
            }
        }
        return near_idx;
        };
    edge_vector cross_edges;                                                    // 交叉边集
    res_points = std::move(loops_points);                                       // 将loop的点集全部移动结果点集中
    size_t loops_edges_siz = loops_edges.size();                                // 记录loop边集的结尾下标
    res_edges = std::move(loops_edges);                                         // 移动到结果边集中
    const int loops_siz = res_points.size();                                    // 记录点集结尾下标
    res_points.insert(res_points.end(), new_points.begin(), new_points.end());  // 将已经筛选过一轮的面上采样点集添加到结果点集
    // 以下lambda用于处理水密性问题
    auto process_result_edges = [&](edge_vector& edges, std::unordered_map<CDT::Edge, std::vector<double>> intersections, auto get_sweep_cord, auto get_intersection_point) {
        for (int i = 0; i < edges.size(); i++) {
            auto [idx_1, idx_2] = edges[i].verts();  // 获取边上的两个点的索引
            if (new_map[idx_1] != NOTFOUND && new_map[idx_2] != NOTFOUND) {
                // 如果两个点都在面内
                if (intersections.find(edges[i]) == intersections.end()) res_edges.push_back(CDT::Edge(new_map[idx_1] + loops_siz, new_map[idx_2] + loops_siz));  // 插入这个边
            }
            else if (new_map[idx_1] != NOTFOUND || new_map[idx_2] != NOTFOUND) {                                                                                // 如果有一个点在面内，另一个点不在面内，说明与loop产生交点
                if (new_map[idx_1] == NOTFOUND) {                                                                                                                 // 将面内的点设置为idx_1，方便统一处理
                    std::swap(idx_1, idx_2);
                }
                auto& cords = intersections[edges[i]];  // 获取这个边与所有loop产生的所有交点
                // 因为边集都是垂直的或者水平的，所以交点只需要用一个维度即可表示。类型为double
                if (cords.size() % 2 == 0) {
                    continue;
                }
                const double sweep_cord = get_sweep_cord(face_points[idx_1]);
                point temp = get_intersection_point(sweep_cord, cords.front());  // 获取交点坐标
                auto near_idx = find_near_idx_from_kdtree(temp);
                if (new_map[idx_1] + loops_siz != near_idx) res_edges.emplace_back(new_map[idx_1] + loops_siz, near_idx);  // 连接新的边
                for (int i = 1; i < cords.size(); i += 2) {
                    auto temp1 = get_intersection_point(sweep_cord, cords[i]);
                    auto near_idx_1 = find_near_idx_from_kdtree(temp1);
                    auto temp2 = get_intersection_point(sweep_cord, cords[i + 1]);
                    auto near_idx_2 = find_near_idx_from_kdtree(temp2);
                    if (near_idx_1 != near_idx_2) res_edges.emplace_back(near_idx_1, near_idx_2);
                    // debug_graph.AddconstraintEdge({res_points[near_idx_1].x, res_points[near_idx_1].y}, {res_points[near_idx_2].x, res_points[near_idx_2].y});
                }
            }
            else {
                auto cords = intersections[edges[i]];
                if (cords.size() % 2) {
                    continue;
                }
                const double sweep_cord = get_sweep_cord(face_points[idx_1]);
                for (int i = 0; i < cords.size(); i += 2) {
                    auto temp1 = get_intersection_point(sweep_cord, cords[i]);
                    auto near_idx_1 = find_near_idx_from_kdtree(temp1);
                    auto temp2 = get_intersection_point(sweep_cord, cords[i + 1]);
                    auto near_idx_2 = find_near_idx_from_kdtree(temp2);
                    if (near_idx_1 != near_idx_2) res_edges.emplace_back(near_idx_1, near_idx_2);
                    // debug_graph.AddconstraintEdge({res_points[near_idx_1].x, res_points[near_idx_1].y}, {res_points[near_idx_2].x, res_points[near_idx_2].y});
                }
            }
        }
        };
    process_result_edges(v_edges, v_intersection_results, GetX, [](const double sweep_cords, const double cords) { return point{ sweep_cords, cords }; });
    process_result_edges(h_edges, h_intersection_results, GetY, [](const double sweep_cords, const double cords) { return point{ cords, sweep_cords }; });

    CDT::RemoveDuplicatesAndRemapEdges(res_points, res_edges);  // 移除重复的边并重新映射边的索引
    for (int i = 0; i < loops_edges_siz; i++) {                  // 遍历标记loop边集,用于后续CDT过滤处理
        markedEdge.insert({ res_edges[i].v1() + 3, res_edges[i].v2() + 3 });
    }
}
void spline_faceter::attachMeshByQuad(QuadTree& tree) {
#ifdef DEBUG
    debug_graph.init(f);
    for (auto& x : tree.points_) {
        debug_graph.AddSamplePoint(std::make_pair(x.x, x.y));
    }
    for (auto& x : points) {
        debug_graph.AddBoundaryPoint(std::make_pair(x.x, x.y));
    }
    for (auto& x : tree.h_edges_) {
        // debug_graph.AddconstraintEdge(std::make_pair(tree.points_[x.v1()].x, tree.points_[x.v1()].y), std::make_pair(tree.points_[x.v2()].x, tree.points_[x.v2()].y));
    }
    for (auto& x : tree.v_edges_) {
        // debug_graph.AddconstraintEdge(std::make_pair(tree.points_[x.v1()].x, tree.points_[x.v1()].y), std::make_pair(tree.points_[x.v2()].x, tree.points_[x.v2()].y));
    }
    for (auto& x : edges) {
        debug_graph.AddconstraintEdge(std::make_pair(points[x.v1()].x, points[x.v1()].y), std::make_pair(points[x.v2()].x, points[x.v2()].y));
    }
    debug_graph.output();
#endif
    ///
    ///
    ///

     // 标记边集，用于标记loop边，用于CDT的限制边过滤
    CDT::EdgeUSet markedEdge;
    edge_vector final_edges;
    point_vector final_points;
    TIMERSTART(ProcessFaceAndLoopsIntersection)

        // 处理面和Loop之间的边界
        ProcessFaceAndLoopsIntersection(tree, points, edges, final_points, final_edges, markedEdge);
    RemoveDuilplicatePointsAndEdges(final_points, final_edges);  // 移除重复的点和边
    TIMEREND(ProcessFaceAndLoopsIntersection)
        // 创建CDT三角剖分对象
        TIMERSTART(CDT_TRI)
        CDT::Triangulation<double>* pcdt = ACIS_NEW CDT::Triangulation<double>;  // CDT进行三角剖分
    // 插入顶点
    pcdt->insertVertices(final_points);
    // 插入边
    pcdt->insertEdges(final_edges);
    auto&& mp = CDT::EdgeToPiecesMapping(pcdt->pieceToOriginals);
    for (auto& x : mp) {
        if (markedEdge.find(x.first) != markedEdge.end()) {
            for (auto& v : x.second) {
                markedEdge.insert(v);
            }
        }
    }
    pcdt->eraseOuterTrianglesAndHoles(markedEdge);
    TIMEREND(CDT_TRI)
        // 对三角形进行后处理
        TIMERSTART(PostProcessTriangle)
        PostProcessTriangle(f, nt, st, pcdt);
    TIMEREND(PostProcessTriangle)

        TIMERSTART(Other)
        // 获取三角形和顶点数据
        const CDT::TriangleVec triangles = std::move(pcdt->triangles);
    const std::vector<CDT::V2d<double>> vertice = std::move(pcdt->vertices);
    // 创建索引网格对象
    const surface& suf = f->geometry()->equation();
    std::vector<UV_VERTEX> uv_vertices;
    for (int i = 0; i < vertice.size(); i++) {
        CDT::V2d<double> pp = vertice[i];
        SPApar_pos param(pp.x, pp.y);
        // 计算顶点位置和法线
        const SPAposition pos = suf.eval_position(param);
        SPAunit_vector normal = suf.eval_normal(param);
        // 如果面的方向相反，则反转法线
        if (f->sense()) normal = -normal;
        // 添加顶点到网格
        uv_vertices.emplace_back(pos, normal, param);
    }
    NewIndexedMesh(f, uv_vertices, triangles);
    ///
    ///
    ///
#ifdef DEBUG
    auto Edges_set = CDT::extractEdgesFromTriangles(triangles);
    debug_graph.AddTriangles(vertice, Edges_set);
    debug_graph.output();
#endif
}

outcome spline_faceter::facet() {
    // 初始化
    TIMERSTART(TOTAL);
    init();
    TIMERSTART(QUADTREE);
    SPAinterval u_interval = f->geometry()->equation().param_range_u();
    SPAinterval v_interval = f->geometry()->equation().param_range_v();
    // 初始化叶子节点的边界框
    Box box = Box{
      {u_interval.start_pt(), v_interval.start_pt()},
      {u_interval.end_pt(),   v_interval.end_pt()  }
    };
    QuadTree tree(f, box, nt, st);
    TIMEREND(QUADTREE);

    auto generator = [&](SPApar_pos begin, SPApar_pos end) {
        double min_x = std::min(begin.u, end.u);
        double max_x = std::max(begin.u, end.u);
        double min_y = std::min(begin.v, end.v);
        double max_y = std::max(begin.v, end.v);
        bool is_u = begin.u == end.u;
        bool is_less = is_u ? (begin.v < end.v) : (begin.u < end.u);
        std::vector<SPApar_pos> tmp;
        for (auto& tmp_point : tree.points_) {
            if (is_u) {
                if (tmp_point.x >= min_x && tmp_point.x <= max_x && tmp_point.y > min_y && tmp_point.y < max_y) {
                    tmp.emplace_back(tmp_point.x, tmp_point.y);
                }
            }
            else {
                if (tmp_point.x > min_x && tmp_point.x < max_x && tmp_point.y >= min_y && tmp_point.y <= max_y) {
                    tmp.emplace_back(tmp_point.x, tmp_point.y);
                }
            }
        }
        if (is_u) {
            sort(tmp.begin(), tmp.end(), [&](const SPApar_pos& a, const SPApar_pos& b) { return is_less ? a.u < b.u : a.u > b.u; });
        }
        else {
            sort(tmp.begin(), tmp.end(), [&](const SPApar_pos& a, const SPApar_pos& b) { return is_less ? a.v < b.v : a.v > b.v; });
            return tmp;
        }
        };

    TIMERSTART(LOOPGAIN);
    auto loop_processor = LoopProcessor(f);
    loop_processor.RegisterGenerator(generator);
    loop_processor.Execute();
    for (auto x : loop_processor.points) {
        points.emplace_back(x.u, x.v);
    }
    for (int loop_idx = 0; loop_idx < loop_processor.loops_indexs.size(); loop_idx++) {
        auto& point_indexs = loop_processor.loops_indexs[loop_idx];
        int point_indexs_size = point_indexs.size();
        if (point_indexs_size < 2) continue;
        for (int i = 1; i < point_indexs_size; i++) {
            edges.emplace_back(point_indexs[i], point_indexs[i - 1]);
        }
        edges.emplace_back(point_indexs.front(), point_indexs.back());
    }
    TIMEREND(LOOPGAIN);

    TIMERSTART(AttachMeshByQuad);
    attachMeshByQuad(tree);
    TIMEREND(AttachMeshByQuad);
    TIMEREND(TOTAL);
    return outcome();
}

outcome facet_face_spline(FACE* f, double nt, double st) {
    spline_faceter faceter(f, nt, st);
    return faceter.facet();
}
void spline_faceter::decideUVlen() {
}
// 获取periphery点集, 其他loop点集
bool spline_faceter::getLoops() {
    return false;
}
